def plot_dist(data,
              sat_var,
              clim,
              extent=[128.95, 129.5, 35.3, 35.75],
              stn_list=None,
              preset='Ulsan',
              figsize=(10, 10),
              grid_fontsize=13):
    import numpy as np
    from air_toolbox import plot, custom_colorbar
    import matplotlib.pyplot as plt
    import proj_util

    plot.set_figure(fontsize=20)

    _cmap = custom_colorbar('GEMS')
    _fig, _ax = plot.map_plot(preset=preset,
                              grid_bool=False,
                              road=True,
                              road_region='Ulsan',
                              extent=extent,
                              figsize=figsize,
                              grid_fontsize=grid_fontsize)

    _im = _ax.pcolormesh(
        data.lond,
        data.latd,
        np.nanmean(data.__dict__[sat_var], axis=2),
        cmap=_cmap,
    )

    if stn_list is not None:
        if type(stn_list) == str:
            stn_list = [stn_list]
        for stn_n in stn_list:
            stn, vaa_list = proj_util.gbrs_info(stn_n, return_vaa=True)
            plot.plot_stn(_ax, stn, c='k', vaa_bool=True, vaa_list=vaa_list)

    _im.set_clim(clim)
    _cbar = plot.colorbar(_fig, _ax, _im, 'horizontal', pad=0.5)

    return _fig, _ax, _cbar


def plot_param(prod):
    if prod == 'NO2':
        plot_prod = 'NO$_2$'
        unit = 'molec. cm$^{-2}$'
    return plot_prod, unit


def label_param(sat_name):
    if sat_name == 'GEMS':
        _sat_colname = 'GEMS_V3'
        _sat_labelname = 'GEMS v3.0'
    elif sat_name == 'GEMS_IUP':
        _sat_colname = sat_name
        _sat_labelname = 'GEMS-IUP'
    elif sat_name == 'combine':
        _sat_colname, _sat_labelname = None, None
    return _sat_colname, _sat_labelname


def plot_scatter(df,
                 sat_var,
                 lim,
                 density=True,
                 option=None,
                 grid=True,
                 prod='NO2',
                 sat_name='GEMS',
                 text=True,
                 inst='Pandora',
                 dpi=20,
                 c='k',
                 tick_fontsize=13,
                 label_fontsize=13):
    from air_toolbox import plot
    plot.set_figure(fontsize=tick_fontsize)

    plot_prod, unit = plot_param(prod)

    if sat_name == 'GEMS':
        _sat_colname = 'GEMS_V3'
        _sat_labelname = 'GEMS v3.0'
    elif sat_name == 'GEMS_IUP':
        _sat_colname = sat_name
        _sat_labelname = 'GEMS-IUP'

    if sat_var == 'VCD':
        # suffix = 'Trop. VCD'
        suffix = 'TrVCD'
    elif sat_var == 'TCD':
        suffix = 'TCD'
    var_name = f'{plot_prod} {suffix}'
    if type(option) == str and option.upper() == 'VAA_AVG':
        df_sat_varname = f'{_sat_colname}_{sat_var}_VAA_AVG'
    else:
        df_sat_varname = f'{_sat_colname}_{sat_var}'
    if density:
        fig, ax = plot.scatter_density(df[f'RS_{sat_var}'],
                                       df[df_sat_varname],
                                       lim,
                                       offsettext=True,
                                       text=text,
                                       dpi=dpi)
    else:
        fig, ax = plot.scatter(df[f'RS_{sat_var}'],
                               df[df_sat_varname],
                               lim=lim,
                               offsettext=True,
                               text=text,
                               tick_fontsize=tick_fontsize,
                               c=c)
    ax.grid(grid)
    ax.set_xlabel(f'{inst} {var_name} [{unit}]', fontsize=label_fontsize)
    ax.set_ylabel(f'{_sat_labelname} {var_name} [{unit}]',
                  fontsize=label_fontsize)
    if label_fontsize > 14:
        offset_text = ax.xaxis.get_offset_text()
        offset_text.set_position((1.2, 0))
    return fig, ax


LENGTH = 10  # 길이를 고정


def calc_coordinate(aa, za):
    """
    azimuth와 zenith angle을 받아 3D 좌표를 계산하여 반환.
    
    Parameters:
        aa (float): azimuth angle (도 단위)
        za (float): zenith angle (도 단위, 0~90도)
    
    Returns:
        x, y, z: 3D 좌표
    """
    import numpy as np
    # 각도를 라디안으로 변환
    _aa_rad, _za_rad = np.radians(aa), np.radians(za)

    # 3D 좌표 계산
    _x = LENGTH * np.sin(_aa_rad) * np.cos(_za_rad)
    _y = LENGTH * np.cos(_aa_rad) * np.cos(_za_rad)
    _z = LENGTH * np.sin(_za_rad)

    return _x, _y, _z


def plot_azimuth_3d(vaa, saa, vza, sza, inst_vaa=None, inst_vza=0):
    """
    위성(vaa, vza)과 태양(saa, sza)의 azimuth & 고도 각도를 받아 3D 시각화하여 출력.
    추가적으로 inst_vaa, inst_vza 값을 받아 초록 점선 화살표로 표시.
    
    Parameters:
        vaa (float): 위성의 azimuth angle (도 단위)
        saa (float): 태양의 azimuth angle (도 단위)
        vza (float): 위성의 고도각 (도 단위, 0~90도)
        sza (float): 태양의 고도각 (도 단위, 0~90도)
        inst_vaa (float): 추가적인 azimuth angle (도 단위)
        inst_vza (float, optional): 추가적인 고도각 (도 단위, 기본값=0)
    
    Returns:
        fig, ax: Matplotlib figure와 axis 객체
    """
    import matplotlib.pyplot as plt
    from air_toolbox import plot
    # 위성, 태양, inst의 3D 좌표 계산
    _sat_x, _sat_y, _sat_z = calc_coordinate(vaa, vza)
    _sun_x, _sun_y, _sun_z = calc_coordinate(saa, sza)

    _lim_h = [-LENGTH - 2, LENGTH + 2]  #horizontal range
    _lim_v = [0, LENGTH * 0.7]  #vertical range
    # 그래프 설정
    plot.set_figure(fontsize=13)
    _fig, _ax = plt.subplots(figsize=(8, 8),
                             subplot_kw={'projection': '3d'},
                             constrained_layout=True)

    # 위성, 태양, inst을 점으로 표시
    _ax.scatter(_sat_x, _sat_y, _sat_z, color='blue', s=50, label='Satellite')
    _ax.scatter(_sun_x, _sun_y, _sun_z, color='red', s=50, label='Sun')
    _ax.scatter(0, 0, 0, color='green', s=50, label='Instrument')

    # 원점에서 위성, 태양, inst까지 선 연결
    _ax.plot([0, _sat_x], [0, _sat_y], [0, _sat_z], 'b-', alpha=0.7)
    _ax.plot([0, _sun_x], [0, _sun_y], [0, _sun_z], 'r-', alpha=0.7)
    if inst_vaa is not None:
        _inst_x, _inst_y, _inst_z = calc_coordinate(inst_vaa, inst_vza)
        _ax.quiver(0,
                   0,
                   0,
                   _inst_x,
                   _inst_y,
                   _inst_z,
                   color='green',
                   linestyle='-',
                   arrow_length_ratio=0.1)
    # 수선의 발 계산 (지면과의 교차점)
    _ax.plot([_sat_x, _sat_x], [_sat_y, _sat_y], [0, _sat_z], 'b:')
    _ax.plot([_sun_x, _sun_x], [_sun_y, _sun_y], [0, _sun_z], 'r:')
    if inst_vza != 0:
        _ax.plot([_inst_x, _inst_x], [_inst_y, _inst_y], [0, _inst_z],
                 'g:',
                 label='Instrument Foot')

    # set the range of x, y, z
    _ax.set_xlim(_lim_h)
    _ax.set_ylim(_lim_h)
    _ax.set_zlim(_lim_v)
    # remove ticks
    _ax.set_xticks([])
    _ax.set_yticks([])
    _ax.set_zticks([])
    # set labels
    _ax.set_xlabel("West (-) / East (+)", fontsize=13)
    _ax.set_ylabel("South (-) / North (+)", fontsize=13)

    # 중앙 십자 표시 (지면에 교차선 추가)
    _ax.plot([-LENGTH, LENGTH], [0, 0], [0, 0], 'k-', alpha=0.7)
    _ax.plot([0, 0], [-LENGTH, LENGTH], [0, 0], 'k-', alpha=0.7)

    # 동서남북 라벨 추가
    _ax.text(0, LENGTH + 1, 0, 'N', fontsize=12, ha='center', va='bottom')
    _ax.text(LENGTH + 1, 0, 0, 'E', fontsize=12, ha='left', va='center')
    _ax.text(0, -LENGTH - 1, 0, 'S', fontsize=12, ha='center', va='top')
    _ax.text(-LENGTH - 1, 0, 0, 'W', fontsize=12, ha='right', va='center')

    _ax.view_init(elev=30, azim=-45)

    # 격자 및 범례
    _ax.legend()
    _ax.grid(True, linestyle='--', alpha=0.6)

    return _fig, _ax


def diurnal_plot(group_obj,
                 sat_var,
                 inst_name='P150',
                 figsize=(15, 5),
                 label_fontsize=20,
                 tick_fontsize=15,
                 legend_fontsize=28,
                 hist=True,
                 threshold=True,
                 wind=False):

    import matplotlib.pyplot as plt
    import matplotlib.dates as mdates
    from air_toolbox.plot import set_figure
    import numpy as np
    set_figure(fontsize=tick_fontsize)
    _fig, _ax = plt.subplots(figsize=figsize, constrained_layout=True)

    _ax.patch.set_visible(False)
    ### plot median
    _ax.errorbar(group_obj.median.index,
                 group_obj.median[f'GEMS_V3_{sat_var.upper()}'],
                 yerr=[
                     group_obj.median[f'GEMS_V3_{sat_var.upper()}'] -
                     group_obj.q1[f'GEMS_V3_{sat_var.upper()}'],
                     group_obj.q3[f'GEMS_V3_{sat_var.upper()}'] -
                     group_obj.median[f'GEMS_V3_{sat_var.upper()}']
                 ],
                 fmt='^r-',
                 label='GEMS',
                 elinewidth=2,
                 capsize=7,
                 capthick=3,
                 ms=15,
                 zorder=700)
    _ax.errorbar(group_obj.median.index,
                 group_obj.median[f'RS_{sat_var.upper()}'],
                 yerr=[
                     group_obj.median[f'RS_{sat_var.upper()}'] -
                     group_obj.q1[f'RS_{sat_var.upper()}'],
                     group_obj.q3[f'RS_{sat_var.upper()}'] -
                     group_obj.median[f'RS_{sat_var.upper()}']
                 ],
                 fmt='-ok',
                 label=inst_name,
                 elinewidth=2,
                 capsize=7,
                 capthick=3,
                 ms=15,
                 zorder=600)
    if wind:
        _ax_wind = _ax.twinx()
        _ax_wind.spines['right'].set_position(('outward', 90))
        _ax_wind.errorbar(group_obj.median.index,
                          group_obj.median.wspd10,
                          yerr=[
                              group_obj.median.wspd10 - group_obj.q1.wspd10,
                              group_obj.q3.wspd10 - group_obj.median.wspd10
                          ],
                          fmt='-sb',
                          label='Wind Speed',
                          elinewidth=2,
                          capsize=7,
                          capthick=3,
                          ms=15,
                          zorder=500)
        _ax_wind.set_ylabel('Wind Speed (m/s)', fontsize=28, color='blue')
        _ax_wind.tick_params(axis='y', labelsize=25, colors='blue')
        _ax_wind.set_ylim([0, 8])

    _ax.set_xlabel('Time [KST]', fontsize=label_fontsize, labelpad=10)
    _ax.set_ylabel(f'NO$_2$ TrVCD [molec. cm$^{{-2}}$]',
                   fontsize=label_fontsize,
                   labelpad=10)
    _ax.set_ylim([0, 2e16])
    _ax.set_xlim([7.5, 17.5])
    _ax.tick_params(axis='x', which='major', length=10, width=2)
    _ax.set_xticks(np.arange(8, 18, 1))

    if hist:
        _ax2 = _ax.twinx()
        _ax.set_zorder(_ax2.get_zorder() + 1)
        ### plot histogram of data count
        _ax2.hist(group_obj.Hour,
                  bins=np.arange(group_obj.count.index.min(),
                                 group_obj.count.index.max() + 2) - 0.5,
                  color='green',
                  alpha=0.3,
                  rwidth=0.5)
        _ax2.set_ylabel('Data count', fontsize=label_fontsize, labelpad=10)
        _ax.legend(loc='upper right', fontsize=legend_fontsize, ncols=2)

    else:
        for _key, _group in group_obj.groupby:
            _ax.text(_key,
                     1.8e16,
                     f'{_group.shape[0]}',
                     fontsize=15,
                     color='green',
                     horizontalalignment='center')
            _ax.legend(fontsize=legend_fontsize,
                       ncols=2,
                       bbox_to_anchor=(1.01, 1.2))
            _ax.grid(True)

    return _fig, _ax


class hourly_group:

    def __init__(self, df, sat_var, sat_name):
        self.colprefix, self.labelname = label_param(sat_name)
        if sat_name == 'combine':
            _colname = [
                f'GEMS_V3_{sat_var.upper()}', f'GEMS_IUP_{sat_var.upper()}',
                f'RS_{sat_var.upper()}', 'Hour'
            ]
        else:
            _colname = [
                f'{self.colprefix}_{sat_var.upper()}', f'RS_{sat_var.upper()}',
                'Hour'
            ]
        _df = df[_colname]
        self.Hour = _df['Hour']
        self.groupby = _df.groupby('Hour')
        self.mean = self.groupby.mean()
        self.std = self.groupby.std()
        self.count = self.groupby.count()
        self.median = self.groupby.median()
        self.max = self.groupby.max()
        self.min = self.groupby.min()
        self.q1 = self.groupby.quantile(0.25)
        self.q3 = self.groupby.quantile(0.75)
